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^\ • A stochastic model of autocatalytic chemical reactions is studied both numerically and analyt- 

£N) ' ically. The van Kampen perturbative scheme is implemented, beyond the second order approxi- 

mation, so to capture the non Gaussianity traits as displayed by the simulations. The method is 
targeted to the characterization of the third moments of the distribution of fluctuations, originating 
from a system of four populations in mutual interaction. The theory predictions agree well with 
^ , the simulations, pointing to the validity of the van Kampen expansion beyond the conventional 

Gaussian solution. 
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^ 1 I. INTRODUCTION 

The cell is a complex structural unit, that defines the building block of living systems [1]. It is made of by a tiny 
membrane, constituted by a lipid bilayer, which encloses a finite volume and protects the genetic material stored 
inside. The membrane is semi-permeable: nutrients can leak in and serve as energy storage to support the machinery 
functioning. Metabolism converts energy into molecules, i.e. building cell components, and releases by-product. 

Evolution certainly guided the ancient supposedly minimalistic cell entities, the so-called protocells 2-4], through 
subsequent steps towards the delicate and complex biological devices that we see nowadays. Focusing on primordial 
cell units, back at the origin of life, the most accredited scenario dictates that chemical reactions occurred inside 
(/*■) | vesicles, small cell-like structures in which the outer membrane takes the form of a lipid bilayer Vesicles possibly 
■ ' ■ defined the scaffold of prototypical cell models, while it is customarily believed that autocatalytic reactions might 
have been at play inside primordial protocell. The shared view is that protocell's volume might have been occupied by 
. interacting families of replicators, organized in autocatalytic cycles. A chemical reaction is called autocatalytic if one 
7-H ' of the reaction products is itself a catalyst for the chemical reaction. Even if only a small amount of the catalyst is 
present, the reaction may start off slowly, but will quickly develop once more catalyst is produced. If the reactant is not 
replaced, the process will again slow down producing the typical sigmoid shape for the concentration of the product. 
All this is for a single chemical reaction, but of greater interest is the case of many chemical reactions, where one or 
more reactions produce a catalyst for some of the other reactions. Then the whole collection of constituents is called an 
autocatalytic set. Autocatalytic reactions have been invoked in the context of studies on the origin of life as a possible 
solution of the famous Eigen's paradox [j| . This is a puzzling logic concept which limits the size of self replicating 
molecules to perhaps a few hundred base pairs. However, almost all life on Earth requires much longer molecules to 
encode their genetic information. This problem is handled in living cells by the presence of enzymes which repair 
mutations, allowing the encoding molecules to reach large enough sizes. In primordial organisms, autocatalytic cycles 
might have contributed to the inherent robustness of the system, translating in a degree of microscopic cooperation 
that successfully prevented the Eigen's evolutionary derive towards self-destruction to occur. It is therefore of interest 
to analyze the coupled dynamics of chemicals organized in extended cycles of autocatalytic reactions. 

It is in this context that our work is positioned. We will in particular consider a model of autocatalytic reactions 
confined within a bounded region of space. The model was pioneered by Togashi and Kaneko 6] and more recently 
revisited by @, H[. It was in particular shown that fluctuations stemming from the intimated discreteness of the 
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scrutinized medium can seed a resonant effect yielding to organized macroscopic patterns, both in time and space 
0- 

As we shall clarify in the forthcoming discussion, the model here examined is intrinsically stochastic and falls in 
the realm of the so called individual-based description. The microscopic dynamics follows explicit rules governing 
the interactions among individuals and with the surrounding environment. Starting from the stochastic scenario 
and performing a perturbative development (van Kampen expansion 0) with respect to a small parameter which 
encodes the amplitude of finite size fluctuations, one obtains, at the leading order, the mean-field equations, i.e. the 
idealized continuum description for the concentration amount. These latter govern in fact the coupled evolution of the 
average population amount, as in the spirit of the deterministic representation. Including the next-to-leading order 
corrections, one obtains a description of the fluctuations, as a set of linear stochastic differential equations. Such a 
system can be hence analyzed exactly, so allowing us to quantify the differences between the stochastic formulation 
and its deterministic analogue. This analysis was performed in Q with reference to the a-spatial version of model, 
and in Q where the notion of space is instead explicitly included. 

In this paper, we take one step forward by analytically characterizing the fluctuations beyond the second order in the 
van Kampen perturbative scheme [jl H3] , i-e. the Gaussian approximation, and so quantifying higher contributions in 
the hierarchy of moments of the associated distribution. As we shall demonstrate, and with reference to the analyzed 
case study, we can successfully quantify non Gaussian fluctuations, within the van Kampen descriptive scenario, in 
agreement with the recent investigations of Grima and collaborators and previous indications of Risken and 
Vollmer [lj]. 

Again, let us emphasize that fluctuations do not arise from an external noise. Despite the evidence that it is always 
present in actual population dynamics and that it is an essential ingredient of life processes, noise is often omitted. 
When instead considered, it is frequently assumed to act as a source of disorder and it is included in the dynamics 
as an external elements. At variance, the individual-level approach allows us to investigate the unavoidable intrinsic 
noise, which originates from the discreteness of the system and that has to be considered in any sensible model of 
natural phenomena. 

The paper is organized as follows. In the following section we will introduce the model under scrutiny. Then we will 
turn to discussing the associated master equation, derive the mean field equation, and characterize the fluctuations 
within the Gaussian approximation. Non Gaussian traits are revealed via numerical (stochastic) simulations for small 
system sizes. These features are analytically inspected and explained in section VI by working in the framework of a 
generalized Fokker-Planck formulation where the role of the finite population is explicitly accommodated for. Finally, 
in section VII we sum up and conclude. 



II. THE MODEL 

The autocatalytic reaction scheme as introduced in [7| describes the dynamics of k species which interact according 
to the following rules 

X-i + Xi + i — > 2X; + i with Afc + i = X\ 
E ^ X t 

X % A E (1) 

where denotes an element of the i-th species, while E is the null constituent or vacancies. The parameter r» (with 
rfe + i = n) is the autocatalytic process rate, while oii and j3i are the rates at which the molecules appear and disappear 
from the system. The size of the system is denoted by N, then rii + tie = N, where tie is the number of E. 

It is worth emphasizing that the concept of vacancies E enables us to accommodate for a finite carrying capacity 
of the hosting volume. The approach can be readily extended to the case where the space is accounted for by 
formally dividing the volume in small patches, each being characterized by a limited capacity. Species can then 
migrate between neighbors cells, therefore visiting different regions of the spatial domain in which they are confined. 
This generalization is discussed in We will here solely consider the a-spatial version of the model, aiming at 
characterizing the fluctuations beyond the canonical Gaussian approximation. In the following, we will introduce the 
master equation that rules the stochastic dynamics of the system defined by the closed set of chemical equations |T]) . 



III. THE MASTER EQUATION AND ITS EXPANSION 



Let us start by introducing the master equation that governs the evolution of the stochastic system described 
above. First, it is necessary to write down the transition rates T(n'\n) from the state n to the state «/, where 
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n = (ni, ...,rik) is the vector whose components define the number of elements of each species at time t. These 
transition rates are 

T(ni, ...,rn - l,n i+ i + 1, . . . ,nk\n) = r i+ i- 




T(ni, ...,rii-l,.. .,n k \n) = fa-^ 



In this way, the differential equation for the probability P(n, t) reads 

d -P(n, t) = ^(e+er +1 )T(ni, ...n* - 1, n l+1 + 1, ...n k )P(n, t)+ 

i—l 

k 

J2( £ i ~ l)T(n x , rn + 1, n i+1 , n k )P(n, t)+ (2) 

i=l 
k 

- l)r(m, Hi - 1, n i+ i, n fc )P(n, t) 



i=l 



where ef are the step operators which act on an arbitrary function /(x) as e^/(x) = /(. . . , Xi ± 1, . . .). 

The above description is exact: no approximations have yet been made. At this stage we could resort to numerical 
simulations of the underlying chemical reactions by means of the Gillespie algorithm |13l . Il4j . This method produces 
realizations of the stochastic dynamics which are formally equivalent to those found from the master equation 
Averaging over many realizations enables us to calculate quantities of interest. We will comment on the results of such 
simulations, in the following. A different route is however possible which consists in drastically simplifying the master 
equation, via a perturbative calculation, the van Kampen system size expansion 0, [To| . It is effectively an expansion 
in powers of A -1 / 2 , which to the leading order (N — > oo) gives the deterministic equations describing the system, 
while at next-to-leading order returns the finite N corrections to these. The method consists in putting forward the 
ansatz: 



a^ + T/F 



(3) 

where is the «-th component of the fc-dimensional stochastic variable £ = (£1,62, •■•)■ To proceed in the analysis 
we make use of the working ansatz (j3]) into the master equation ([2]). Then, it is straightforward to show that the 
operator ef can be approximated as: 

± . , 19 Id 2 Id 3 
ef = 1 ± + — ± 



AV2^ 2A<9£ 2 3!A 3 / 2 9£3 

The first step in the perturbative calculation consists in expliciting in the master equation the dependence on the 
concentration vector y — n/A. It is legitimate to assume that this latter quantity changes continuously with time, 
as far as each instantaneous variation is small when compared to the system size. We therefore proceed by defining 
the following distribution: 



II($,i)=P(y,i)=P(V(i) + ^,i). 

A simple manipulation yields to: 

op ^-Aan#i an 

= -VA> — A . 

dt ^ d& dt dt 

Similarly one can act on the right hand side of Eq. ([2]) and hierarchically organize the resulting terms with respect 
to their A-dependence. The outcome of such algebraic calculation are reported in the following. We will in particular 
limit our discussion to the Gaussian approximation, by neglecting, at this stage, the N~ 3 / 2 terms. We will then return 
on this important issue and discuss the specific role that is played by iV~ 3 / 2 corrections. 
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FIG. 1: Temporal evolution of one of the species concentrations for a system composed by 4 species and parameters set as 
N = 8190, n = 10 and on — Pi = 1/64 Vi. The noisy (red online) line represents one stochastic realization thought the Gillespie 
algorithm [13, EH, while the dashed black line shows the numerical solution of the deterministic system given by Eq. ((4]). 



A. The N~ 1/2 terms 
As concerns the terms of order N~ i one obtains: 

k k k k 

i v-^ du d(j>i i v-^ , i , i \ <9n i r n i ( sr^ . ^ i 



where the rescaled time r is defined as t = t/N. Thus the following system of differential equations holds for the 
concentration amount fa 

-^r = nfa-ifa - n+ifafa+i +«! ^ m I ~ ^ h ( 4 ) 

which in turn corresponds to working within the so-called mean field approximation and eventually disregard finite 
size corrections. We should emphasize that Eqs. ((4]) are obtained by elaborating on the exact stochastic chemical 
model and exploring the limit for infinite system size N — > oo. 

To make contact with previous investigations Q we shall assume the simplifying setting with /3, = ft, a, = a and 
Ti = r Vi. Under this condition, all species asymptotically converge to the fixed point fa which is readily calculated 
as: 




a(l 0* = ^. (5) 



m — 1 



We now turn to numerical simulation based on the Gillespie algorithm and discuss the case with k = 4 species. As 
reported in Fig. [TJ once the initial transient has died out, the numerically recorded time series keep on oscillating 
around the reference value as specified by relation §5§ . The mean field dynamics has conversely relaxed to the deputed 
equilibrium value. These oscillations stem from the finite size corrections to the idealized mean field dynamics and 
will be inspected in the following. We will be in particular concerned with characterizing the statistical properties of 
the observed signal, and quantify via rigorous analytical means the moments of the distribution of the fluctuations. 



B. The N 1 corrections 

Finite size effects related to the N^ 1 corrections result in the Fokker-Planck equation: 

1=1 3=1 1=1 J 
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which governs the evolution of the distribution II(-). Here Aj(£) reads: 



r i+10i+l)£i 



k 



£m ft£z 



while bij stands for the element ij of matrix B defined as: 




-1 +Oi(l- E m =l 0m) + ft 



if i-- 

if i = 

if 3 - 

if |i 



= 4+1 
: i - 1 



For the sake of clarity we shall introduce the matrix M of elements m l3 



and so rewrite A, as: 



—Oil 



defined as: 



if i = j 

if j = i + 1 

if j = i — 1 

if |t-i|>l 



i' = l 

The Fokker-Planck equation ([6]) has been previously obtained in Q and shown to explain the regular oscillations 
displayed in direct stochastic simulations. The oscillations, in fact, materialize in a peak in the power spectrum of 
fluctuations which can be analytically calculated working in the equivalent context of the Langevin equation. Here, 
we take a different route and reconstruct the distribution of fluctuations through the calculation of the associated 
moments. To allow for analytical progress, we will assume again identical chemical reactions rates for all species, 
namely r; = r, /3j = (3 and a; = a Vi. Moreover, we will focus on the fluctuations around the equilibrium and so 
require </>,; = 4>* Vi. Under these conditions the matrix M is circulating and can be cast in the form: 





' m 


m i 


TO 2 


TO 2 . 


• m 3 




m 3 




m i 


TO 2 . 


■ m 2 


M = 


m 2 


m 3 


to 


TOl . 


■ m 2 




_ mi 


TO 2 


to 2 


TO 2 . 


■ m 



with toq = —a — ft, m\ 



a, to 2 = —a, and = r<j>* — a. The k x k matrix reads instead: 





r b Q 


h 


.. 


. 


h 




bi 




h .. 


. 





B = 





h 


b 


. 







. h 





.. 


. h 


&o 



where b = 2r<p*<p* + a(l - k(f>*) + (3<p* , and b x = -r^*0*. 

We recall that the solution of the Fokker-Planck equation ^ is a multivariate Gaussian which is univocally 
characterized by the associated families of first and second moments. Working within this setting, it is hence sufficient 
to derive the analytical equations that control the time evolution of the first two moments of the distribution. We will 
in particular provide closed analytical expressions for the asymptotic moments and draw a direct comparison with 
the numerical experiments. 



IV. ANALYTICAL ESTIMATES OF THE FLUCTUATIONS DISTRIBUTION MOMENTS 



Define the moment of order p for £j the quantity 
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Let us illustrate the analytical procedure that is here adopted, with reference to To this end we start from 

Eq. © and multiply it on both sides by the factor £ 2 . Integrating over M. k in d£ — (d£i, <i£ 2 , •••) <^6): yields: 

/ g-^m,r)dt= J g y ££-Mm&T)d£+ [ ^l^^-hjmi^r)^. (7) 

* i 1 ' i.j ■* 

Consider the right hand side of Eq. (O and operate two successive integrations by parts. Just two terms survive, 
as it can be trivially argued for. Hence, bringing out the time derivative from the integral at the left hand side of Eq. 
0, the sought equation for the second moments reads: 

(£?) = 2m M (£?) + 2ra M _ 1 + 2to m+1 (&& +1 ) + 2m. M+2 (£ i & +2 ) + 6 M (8) 

where i = 1, 4. Use has been made of the definitions of the coefficients m^. With analogous steps one immediately 
obtains the differential equation that governs the time evolution of quantity (66) : 

(66) ="v(£i£ i+ i) + m J+M (£ 2 ) +to m+ i(^ 2 +1 ) +tom +2 (&+i&+2> 

+ mi + i,i+i(66+i) + m i+ i )i+2 (&6+2) + ( 9 ) 

+ m i, 2+3(^1+3^+1) + &i,i+l- 

which, in practice, encodes the degree of temporal correlation between species i and j. The picture is completed by 
providing the equations for the first moments which read: 

(6) = m lti {£i) + m M „i(&_i) + (&)m M+ i - m M+2 (c; l+2 ). 

Taking into account all possible permutations of the involved indexes both ranging in the interval from 1 to 4, 
and recalling the Eq.s ©-([9]), one eventually obtains a closed system of ten coupled ordinary differential equations. 
For the simplified case ri = r, oti = a, fa = /3 Vi, this latter can be cast in a compact form by introducing the matrix: 



K = 



By further defining: 

x = [«?) 

and the vector 



/ 2toq 2toi 


2m 2 


2m 3 

















\ 


m 3 


2mo 


mi 


m 2 


mi 


mi 


m 2 


m 3 








m 2 


m 3 


2 mo 


mi 





mi 





m 2 


m 3 





mi 


m 2 


m 3 


2 mo 








mi 





m 2 


m 3 





2m 3 








2mo 


2mi 


2m 2 














TO 2 


m 3 





m 3 


2mo 


mi 


mi 


m 2 








TOl 





m 3 


m 2 


m 3 


2mo 





mi 


m 2 








2m 2 








2m 3 





2mo 


2mi 











mi 


m 2 





m 2 


m 3 


m 3 


2mo 


mi 


V 








2mi 








2m 2 





2m 3 


2m / 


(66) 


(66 


> (66) 


(6 2 > 


(6e 


3) (66) 


(6 2 > 


(66) 


D 


= [60 


hi 


h 


60 


bi 


b 


61 


60] 





(6 2 )] 



one gets 

X = KX + D . 

As anticipated, we focus in particular on the late time evolution of the system, i.e. when the fluctuations' distribution 
has converged to its asymptotic form. This request translates into the mathematical condition X = 0, which implies 
dealing with an algebraic system of equations. Given the peculiar structure of the problem, and by invoking a 
straightforward argument of symmetry |l7| . one can identify three families of independent unknowns, namely: 



(6 2 ) = (6 2 > - (6 2 > - (6 2 > =: Ti 
(66) = (66) = (66) = (66) == r 2 
(66) = (66) =: r 3 . 



(10) 
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FIG. 2: Plots of the moments Fi and F2 as functions of a. The black lines show the theoretical predictions given by Eq. 
while the (colored online) symbols represent the numerical simulations of the stochastic problem. Each symbol corresponds to 
a different component of the family according to (|10[1 . Parameters are set as N = 2000, a = ft. 



Closed analytical expressions for the unknowns Y\, T2 and F3 as a function of the chemical parameters can be derived 
and take the form: 



2&o _ 


h 


5a 


5a 


60 


36i 


10a 


+ 10a 


bo 


_ h 


10a 


5a 



(11) 



In deriving the above, we have assumed a further simplifying condition, namely a = ft. The adequacy of the predictions 
is tested in Fig. [21 where Ti and T2 arc plotted versus the independent parameter a. Recalling the explicit forms of bo 
and 61 one can immediately appreciate that T3 is indeed independent of a. For this reason we here avoid to include T3 
in Fig. f2J One can moreover make use of the knowledge of the moments to reconstruct the profile of the distribution 
II(£). In particular, and due to the symmetry of the model, we solely focus on the marginal distribution IT(£) = II(£j) 
for i = 1,...,4. In practice, we project the distribution in a one-dimensional subspace by integrating over three out 
of four scalar independent variables In Fig. |3l a comparison between theory and stochastic simulations (relative 
to small N values) is drawn. While the agreement is certainly satisfying, deviations from the predicted Gaussian 
profile manifest as the population size shrinks. As we shall demonstrate, these distortions, which materialize in a 
skewed distribution, can be successfully explained within an extended interpretative framework that moves from the 
van Kampen system size expansion. In the following section we will hence extend the calculation beyond the Gaussian 
approximation. In doing so we will operate in the general setting for a 7^ ft, but then specialize on the choice a = ft 
to drastically reduce the complexity of the inspected problem. 



V. BEYOND THE GAUSSIAN APPROXIMATION 



We shall here go back to discussing the higher orders, N~ 3 / 2 corrections to the Fokker-Planck equation. We will 
in particular consider the various terms that contribute to the generalized Fokker-Planck equation grouping them as 
a function of the order of the derivative involved. 
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FIG. 3: Comparison between the stationary marginal Gaussian distribution and the stochastic simulations (the y-axis has a 
logarithmic scale). The solid (red on line) line shows the theoretical prediction according to the van Kampen theory. The 
(green online) circles represent the numerical distribution for a system with N = 200, while the (red online) triangles refer to 
a system with N — 2000. For all the curves r = 10, a = /3 = 0.1. 



The order N 3 / 2 terms that involve the first derivatives can be expressed as: 



n 



d d 
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1&&+iII(£,t) = ^2qT [ r i+lZiZi+l - »"i6-l&]n(£ , t) 
i=l 



k k 



i=l 3 = 1 



where lij are the elements of the k x k circulant matrix L which, for r% — rVi reads: 



L 



( -r Q ... r \ 

r -r 

r 

-r 

r -r 

V -r ... r / 



The order iV 3 / 2 contribution which depends on the second derivatives can be also expressed in a matricial form. 
In fact, we have: 



1 



d 



d 



d 



o 2 



d d 



k k 



i — 1 j — 1 



where the k x k matrix D of elements dij 1 for Ti = r Vi, reads: 



-r(60* + &+i0*) if 5=* + l 

if J =i-l 

otherwise. 



Finally, the third order derivatives contribute as: 



1 d 3 



3! 



a 9 



1 d 3 



d 3 
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where, we introduced the matrix E defined as: 



f r(P) 2 




3 = i 




w 


= i + 


1 






3 = i 


- 1 


w 


= i — 


1 


rm 2 




3 = i 


+ 1 


w 


= i 




-r((j>* 


2 


3 = i 


+ 1 


w 


= i + 


1 


-r{<j>* 


2 


3 = i 




w 


= i — 


1 


-r(cj)* 


2 


3 = i 


- 1 


w 


= i 




{ -w 


-a(l-fc<£*)] 


i = 3 


= w 









In conclusion, one gets the following equation for the distribution II(£, r) [15| : 

1 k d 2 1 k d 3 

+ E jjF-jr' 1 ^^ - mm . E ^7 '-^- < 12 ) 

where: 

He 
3=1 

We will refer to the latter as to the generalized Fokker-Planck equation. In the following section we will discuss the 
corrections to the Gaussian approximation as deduced by the above mathematical framework. 



VI. NON GAUSSIAN CORRECTIONS TO THE MOMENTS OF THE DISTRIBUTION 

Starting from Eq. (fT2|) we shall now assume k = 4 and calculate the first three moments of the asymptotic 
distribution of the fluctuations around the mean field equilibrium. Clearly, the derivation can be in principle extended 
to evaluate the contribution of higher moments. The algebraic complexity of such an extension is however considerable 
and for this reason the third is the largest moment here characterized. The conclusions are nevertheless rather 
interesting as evaluating the third moment allows us to quantify the observed degree of skewness in the distribution 
of fluctuations. 

When it comes to the first moment one gets: 

^(Ci) = mAZi) + n^,i-l(&-l) + + ra M+2 (& +2 ) + -^72pi,i-l(6£i-l) + ii,i+l(&&+l>] (13) 

This equation differs from the one obtained in section llVl for the additional contribution 

Ml/2 ['*,*-! + h,i+l{£i€i+l}] ■ 

Thanks to the symmetry of the system, which ultimately stems from having assumed rj = r Vi, we can operate in a 
highly simplified framework. We notice in fact that the above term is function of the second moments, which have 
been estimated above and quantified as = (CiCi-i) = ^2 + o(l/>/N). Further we observe that l{ t i-x — — 

Hence the corrections to the Gaussian solution as exemplified in Eq. (fl"3|) contribute with an overall term of order 
AT- 3 / 2 , which can be legitimately neglected at this level of approximation. In conclusion the equation for the first 
moments is identical to that obtained in section HVl 

Working in complete analogy, for the second moments we find: 

(e 2 ) = 2m w <6 a ) + 2mi, i+2 {te i+2 ) + 2m ifi+ i<£<6+i) + 2m M _i<&6_i) + b iti + + i^-i^fi-i)] 

for the variance of each involved species (recalling that the first moments are indeed null) and 

= wii,»(£i&+i) + ra M+2 (&+i&+2) + m i>i+1 {(% +1 ) + m i+M (£ 2 ) + m i+M+1 (&£ i+ i) 
+m i+lti+2 (^i +2 ) + m iyi+1 (£f) + m it i +3 (&+1&+3) + 2 bi > i+1 + 2 bi+1,i 
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for the mutual correlation between distinct populations. The index i ranges from 1 to 4. Again the extra contributions 
are limited to the terms stored in square brackets and prove to be negligible at this level of approximation. In fact 
the third order correlations therein involved should scale as TV -1 / 2 as requested by a simple consistency argument 
and as we shall prove a posteriori. Then, also in this case, thanks to the specific form of the matrix L, the additional 
contribution, stemming from third order moments, vanishes. We come hence to the conclusion that the second 
moments are identical to those calculated in the preceding section ITVl working within the Gaussian ansatz. 
Let us now turn to calculating the third moments. After a lengthy derivation we end up with: 

= 3m (a + 3m 3 <£f&-i) + 3m 2 (^ +2 ) + Sm^f^) + 360(6) + ^K(6 2 ) 

+ matefc+i) + ma&fc-i) + 3m 2 <66+2>] + + J—^^) _ 3 r^ +1 )] 

where 777,0 = —2a, mi = —a — r/5, 777 2 = —a, 7773 = —a + r/5, m4 = 2r/5, 7775 = mg = 0, 7777 = — 2r/5, ms — r/25 
and 7779 = —r/5. Here again, and as anticipated in the preceding discussion, we have chosen the simplifying setting 
with a = ft, which consequently implies <jf — 1/5. Elaborating on the symmetry one can identify five families of 
independent moments, which obey to the above and the following differential equations: 

d = 3mo& 2 6-0 + W&ef-i) + 2m,(66-i6 1) + 27773(66-16+2) + ^3(^6+2) 



dt 



+771! (ef) + 777 2 (e 2 6+l) + &0<&-l) + 26l<^i> 
1 

A 

777 8 1 



+-7^72-^3(^1) + "lefe^+i) + m 3 (6+i6-i) + m 2 <6-i6+2) + m 7 (^f)] 



N l/2 t N i/2 



[-2r(6-i£f6+i) + Mti-itt) + m 6-16+2) - r(&£i-i 



and 



-(6 2 6+i) = 3m (6 2 6+i) + 2777 3 (66-i6+i) + 2tt7 1 (66 2 +i) + 2m 2 (66+16+2) 
+777 3 (6 3 ) + mi(6 2 6+ 2 ) + rmitfo-i) + 60(6+1) + 26i(6) 
+ ]Vi72 I TO 3(6 2 +i) + "ie(66+i) + "13(6+16+3) + "12(6+16+2) + «-7 7 (6 2 )] 

+^172 + ^172 + 2r(6-i6+i6 2 ) + r(^+i) - r(6 2 6+i6+2)] 

for adjacent populations with respect to the assumed ordering. For next-to-neighbors correlation one gets: 

^(6 2 e t+2 ) = 3t77 (6 2 6+2) + 2777 3 (66-l6+ 2 ) + 2777l(66+l6+2) + 2T77 2 (66 2 +2) + ™3<6 ? &+l) 

+m!(6 2 6-i> + ™ 2 (6 3 ) + M6+2) + ^72-h^(6 2 ^+i6+2) + r(6-i6+ 2 6 2 >] 

+ Jj~T/2 I TO 2(6 2 +2) + "^4(66+2) + "13(6-16+2) + ^3(6+16+2)] ■ 

Finally, for correlations that involve three distinct species, we find: 

^(66+16-1) = 3t7i (66+i6-i) +w 3 (6+i6 2 -i) + m i(^ 2 +i^-i) +^2(6+26+16-1) 

+ m 3(6 2 6-i) + "11(66+26-1) + "13(66+16+2) + "i2(66 2 -i) + m i(6 2 ^+i) 

+"12(66+1) +M6-1) + M6+1) + ["1-9 (6-16) +"19(6+16) 

+7777(6-16+1)] + ^i72-[-i-(ef+i66-i) +^(6+i66 2 -i)]- 

Clearly, the fourth moments enter the equation for the third ones. To close the system and so enable for quantitative 
predictions, we can estimate the zero-th order contribution to the fourth moments by recalling the Gaussian solution 
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as obtained in Section HVl and neglecting the l/\/~N terms. In formula: 

<£■> = 3((f?)) 2 

The above quantities can be analytically estimated at equilibrium and expressed as a function of respectively I\, T2, 
r 3 , as derived in section PVTl 



A. The asymptotic evolution of the third moments 

Let us now write down the system of differential equations that controls the dynamics of the five independent 
families of moments of order three jl8| . Such a system takes the form 



X = vx + s 

where X is 

x = [($) (^6+2) m+i^i)] 

and the matrix of coefficients V reads: 



V 



Finally the vector S is: 



where: 



(14) 



(15) 



/ 3m 


3toi 


3m 3 


3to 2 









\ 


to 3 


3too 


2TOl + 771 2 


TOl 




2to 3 + 


2to 2 




TOl 


2to 3 + TO 2 


3too 


TO 3 




2toi + 


2to 2 




TO 2 


m 3 


mi 


3too + 2to-2 




2to 3 + 


2 mi 


- m 2 J 


V 


TO 2 + TOl 


TO 3 + TO 2 


TO 3 + TOl 


3too 


+ TOl 4 


- m 3 4 




s = 


= 1/Vn[si 


32 S3 S 4 











si 


= 3to4Fi 


4- 6TO 3 r 2 


4- 3m 2 r 3 






S2 


= "i 3 ri 4 


- TO7T1 4- 


TO 3 r 3 4- TO 2 r 2 - 


f TO 6 r 2 4~ TOg h 


- [-2r(ri) 2 4- 2rr 3 r! - 2rr 3 r 2 4- 2rrir 2 


S3 


= to 3 Ti 4 


- TO 7 ri 4- 


TO 3 r 3 4- TO 2 r 2 - 


- TO 6 r 2 TOg - 


- [2r(ri) 2 - 2rrir 3 + 2rr 3 r 2 - 2rrir 2 ] 


Si 


= TO 4 r 3 4 


- 2m 3 r 2 - 


f TO 2 ri 






S5 


= TO 7 r 3 4 


- 2TO g r 2 









We now turn to numerical simulations to validate the correctness of the theory. Stochastic simulations are performed 
for small systems (N = 200) and the time evolution of the third moments is monitored for each of the considered 
species and by varying the parameter a, while keeping r unchanged. Results are displayed in Fig.s HHH where 
the simulations outcome (symbols) are compared to the theory predictions. The agreement has to be considered 
satisfactory, a conclusion which a posteriori validates the theory assumptions and in particular confirms the predictive 
ability of the van Kampen expansion beyond the Gaussian approximation [ill ] . 
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VII. CONCLUSION 

The study of an extended set of autocatalytic reactions proves interesting in many respects. The system self- 
organizes at the macroscopic level, both in space and time, as follows a non linear resonance mechanism that enhances 
the stochastic fluctuations stemming from the finite size. The spontaneous emergence of collective patterns, as well 
as regular time oscillations in such a system, was recently addressed @, Q by analyzing in detail the underlying 
stochastic process via the celebrated van Kampen expansion, truncated at the Gaussian oder of approximation. In 
Q, it was also speculated that the intrinsic ability of the autocatalytic systems to drive self-organized structures 
might have played a role in the evolutionary selection of efficient cells, starting from minimalistic protocells entities. 
It was in fact argued that oscillatory, spatially extended patterns, might have resonate with the innate ability of a 
vesicle container to divide in two. One could imagine that the oscillations trigger the splitting event and thus favor 
a natural synchronization between the fission of the vesicle and the rate of production of the genetic material stored 
inside, which needs to be passed to the next generation offspring. 

Besides these highly speculative considerations, which deserve to be carefully checked within a self-consistent picture, 
we are here interested in extending the perturbative calculation beyond the second order approximation and challenge 
its adequacy in capturing the deviation from the idealized Gaussian behavior. Recent support on the validity of the 
van Kampen higher orders calculation have been provided by Grima and collaborators [ill ] . We here bring one more 
evidence on the accuracy of the procedure within a rather complex model, where different species are simultaneously 
made to interact. Numerical simulations performed in a stochastic setting with modest sizes of the population 
involved, so to magnify the role played by finite size corrections, confirm the correctness of the theory predictions. 
Due to the complexity of the proposed model, it is not possible to evaluate a large gallery of successive moments 
and so reconstruct the full distribution of fluctuations. The analysis is hence limited to the third moment, which 
however quantifies the degree of skewness of the recorded fluctuations. In a separate contribution 16j, we will return 
on the issue of the validity of the van Kampen ansatz, working within a considerably simpler model that enables 
us to explicitly calculate all the moments of the distribution at any order of the expansion. We are hence able to 
recover a general and exact analytical solution that, we anticipate, agrees very well with the simulations, inline with 
the conclusion of this work. 
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FIG. 4: Plots of Xi (see Eq. I|15j> ) a s functions of the parameter a, for a system with /3 = a, N = 200 and r = 10. The solid 
black lines represent the numerical solution of the system (|14|) . while the symbols refer to the stochastic simulations (each of 
the four symbols is associated to a different species). 




a 

FIG. 5: Plots of X2 as functions of the parameter a. For the parameters' setting and the explanation of the symbols see caption 
of Fig. H 
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FIG. 6: Plots of Xi as functions of the parameter a. For the parameters' setting and the explanation of the symbols see caption 
of Fig. [2 
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FIG. 7: Plots of Xi as functions of the parameter a. For the parameters' setting and the explanation of the symbols see caption 
of Fig. H 
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FIG. 8: Plots of X$ as functions of the parameter a. For the parameters' setting and the explanation of the symbols see caption 
of Fig. a 



